#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "IDoubleMatrix.h"
#include "LatoolException.h"

namespace LatoolNet {
	using namespace System;

	ref class DoubleGeneralMatrix : public IDoubleMatrix {
	internal:
		double * dat;
	private:
		long * ipiv;
		int m_rownum;
		int m_colnum;
		bool m_isFactorized;

		DoubleGeneralMatrix(DoubleGeneralMatrix ^ orig) {
			m_rownum = orig->m_rownum;
			m_colnum = orig->m_colnum;

			dat = new double [orig->m_rownum * orig->m_colnum];
			for (int i = 0; i < m_rownum * m_colnum; i++) {
				dat[i] = orig->dat[i];
			}

			ipiv = new long[m_rownum];
			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}

			m_isFactorized = orig->m_isFactorized;
		};
	public:
		DoubleGeneralMatrix(int rownum, int colnum) {
			m_rownum = rownum;
			m_colnum = colnum;
			
			dat = new double[rownum * colnum];
			for (int i = 0; i < rownum * colnum; i++) {
				dat[i] = 0.0;
			}

			ipiv = new long[rownum];
			
			m_isFactorized = false;
		};

		~DoubleGeneralMatrix(){
			delete[] dat;
			delete[] ipiv;
		};

		virtual IDoubleMatrix ^ Clone() {		
			return (IDoubleMatrix^) gcnew DoubleGeneralMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property double default[int, int] {
			double get(int i, int j) {
				return dat[j * m_rownum + i];
			}
			void set(int i, int j, double value) {
				dat[j * m_rownum + i] = value;
				m_isFactorized = false;
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::DoubleGeneral; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
		};

		virtual void Invert() {

			long n = m_rownum;
			long lda = n;
			long lwork = n;

			double * work = new double[lwork];
			long info = 0;

			dgetri_(&n, dat, &lda, ipiv, work, &lwork, &info);

			delete[] work;

			m_isFactorized = false;

		};

		virtual void Factorize() {

			long m = m_rownum;
			long n = m_colnum;
			long lda = m;
			long info;

			dgetrf_(&m, &n, dat, &lda, ipiv, &info);

			m_isFactorized = true;
		};

		virtual void Solve(IMatrix ^ b) {

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			char trans = 'N';
			long n = m_rownum;
			long nrhs = gb->m_colnum;
			long lda = n;
			long ldb = n;
			long info = 0;

			dgetrs_(&trans, &n, &nrhs, dat, &lda, ipiv, gb->dat, &lda, &info);
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			long n = m_rownum;
			long nrhs = gb->m_colnum;
			long lda = n;
			long ldb = n;
			long info = 0;

			dgesv_(&n, &nrhs, dat, &lda, ipiv, gb->dat, &ldb, &info);

			m_isFactorized = true;

		};

		void SolveLeastSquare(char trans, IMatrix ^ b) {
	
			DoubleGeneralMatrix ^ gb = (DoubleGeneralMatrix^) b;

			long m = m_rownum;
			long n = m_colnum;
			long nrhs = gb->m_colnum;
			long lda = Math::Max(1, m);
			long ldb = Math::Max(1, Math::Max(m, n));
			long lwork = Math::Min(m, n) + Math::Max(Math::Max(1, m), Math::Max(n, nrhs));
			double * work = new double[lwork];
			long info;

			dgels_(&trans, &m, &n, &nrhs, dat, &lda, gb->dat, &ldb, work, &lwork, &info);

			delete[] work;

			if (info == 0) {
				return;
			} else if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			} else {
				throw gcnew LatoolException(String::Format("the {0}-th diagonal element of the triangular factor of A is zero, so that A does not have full rank; the least squares solution could not be computed.", info));
			}
		};

		void SolveSVD(IDoubleMatrix ^ smat, IDoubleMatrix ^ umat, IDoubleMatrix ^ vtmat) {

			char jobu = 'A';
			char jobvt = 'A';
			long m = m_rownum;
			long n = m_colnum;
			long lwork = Math::Max(3 * Math::Min(m, n) + Math::Max(m, n), 5 * Math::Min(m, n));
			double * work = new double[lwork];
			long lda = Math::Max(1, m);
			long ldu = m;
			long ldvt = n;
			long info;

			double * s = new double[Math::Max(1, Math::Min(m_rownum, m_colnum))];
			DoubleGeneralMatrix ^ u = (DoubleGeneralMatrix^) umat;
			DoubleGeneralMatrix ^ vt = (DoubleGeneralMatrix^) vtmat;

			dgesvd_(&jobu, &jobvt, &m, &n, dat, &lda, s, u->dat, &ldu, vt->dat, &ldvt, work, &lwork, &info);

			for (int i = 0; i < Math::Min(m_rownum, m_colnum); i++) {
				smat[i, i] = s[i];
			}
			delete[] work;
			delete[] s;
			
			if (info == 0) {
				return;
			} else if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			} else {
				throw gcnew LatoolException("bidiagonal form B did not converge to zero.");
			}

		}
	};
}